load("DP-repl-data.RData")
library(geepack)
source("IMCI.R")

nboot <- 1000

demthresh <- 0:10
res.post <- res51.post <- as.list(rep(NA, length(demthresh)))
res.post.pearl <- res51.post.pearl <- rep(NA, length(demthresh))
res.post.exog <- res51.post.exog <- matrix(NA, ncol=2, nrow=length(demthresh))
CI.post.pearl <- CI51.post.pearl <- matrix(NA, ncol=2, nrow=length(demthresh))
CI.post.exog <- CI51.post.exog <- matrix(NA, ncol=2, nrow=length(demthresh))

# Estimates for 1886--1992
for(j in 1:length(demthresh)){
    # Create Variables
    D$auth <- ifelse(D$smldmat < demthresh[j], 1, 0)

    # Get bootstrap draws
    res <- matrix(NA, nrow = nboot, ncol = 3, 
                    dimnames = list(1:nboot, c("pearl", "exog.lo", "exog.up")))
    for(b in 1:nboot){
        dds <- unique(D[,"dyad"])
        dyads.b <- sample(dds, length(dds), replace = TRUE)
        index <- c()
        for(k in 1:length(dyads.b)){
            index <- c(index, rownames(D)[D$dyad==dyads.b[k]])
        }
        X <- D[index, "auth"]
        Y <- D[index, "dispute"]
        Z <- D[index, c("smldep", "lcaprat2", "allies", "contigkb", "logdstab", "majpower")]
        dyad <- D[index, "dyad"]
        d <- na.omit(data.frame(Y,X,Z,dyad))
        if (length(unique(X))==1){
            res[b,1] <- NA; res[b,2:3] <- c(NA,NA)
            warning("All obs assigned to one treatment condition at iteration ", b)
        } else {
            fit <- geeglm(Y ~ X + smldep + lcaprat2 + allies + contigkb + logdstab + majpower, 
                family=binomial(link="logit"), data=d, id=dyad, corstr="ar1")
            nd1 <- nd0 <- data.frame(model.matrix(fit)[,-1])
            nd1[,1] <- 1
            nd0[,1] <- 0
            P1 <- predict(fit, nd1, type="response")
            P0 <- predict(fit, nd0, type="response")
            res[b,1] <- mean(pmax(0, P1-P0)/P1, na.rm=T)
            res[b,2:3] <- c(mean(pmax(0, P1-P0)/P1, na.rm=T), 
                mean(pmin(P1, 1-P0)/P1, na.rm=T))
        }
    cat("Threshold ", j, " Resample ", b, " Done. \n")
    }

    res.post[[j]] <- na.omit(res)

    # Point estimates
    fit.post <- geeglm(dispute ~ auth + smldep + lcaprat2 + allies + 
        contigkb + logdstab + majpower, family=binomial(link="logit"), data=D,
        id=dyad, corstr="ar1")
    nd1.post <- nd0.post <- data.frame(model.matrix(fit.post)[,-1])
    nd1.post[,1] <- 1
    nd0.post[,1] <- 0
    P1.post <- predict(fit.post, nd1.post, type="response")
    P0.post <- predict(fit.post, nd0.post, type="response")

    res.post.pearl[j] <- mean(pmax(0, P1.post-P0.post)/P1.post, na.rm=T)

    res.post.exog[j,] <- c(mean(pmax(0, P1.post-P0.post)/P1.post, na.rm=T), 
        mean(pmin(P1.post, 1-P0.post)/P1.post, na.rm=T))

    # Confidence intervals
    CI.post.pearl[j,] <- c(res.post.pearl[j] - qnorm(.975)*sd(res.post[[j]][,1]),
        res.post.pearl[j] + qnorm(.975)*sd(res.post[[j]][,1]))

    CI.post.exog[j,] <- IMCI(l=res.post.exog[j,1], u=res.post.exog[j,2],
        var.l=var(res.post[[j]][,2]), var.u=var(res.post[[j]][,3]))$ci
        
}

# Estimates for 1951--1992
for(j in 1:length(demthresh)){
    # Create Variables
    D51$auth <- ifelse(D51$smldmat < demthresh[j], 1, 0)

    # Get bootstrap draws
    res <- matrix(NA, nrow = nboot, ncol = 3, 
                    dimnames = list(1:nboot, c("pearl", "exog.lo", "exog.up")))
    for(b in 1:nboot){
        dds <- unique(D51[,"dyad"])
        dyads.b <- sample(dds, length(dds), replace = TRUE)
        index <- c()
        for(k in 1:length(dyads.b)){
            index <- c(index, rownames(D51)[D51$dyad==dyads.b[k]])
        }
        X <- D51[index, "auth"]
        Y <- D51[index, "dispute"]
        Z <- D51[index, c("smldep", "lcaprat2", "allies", "contigkb", "logdstab", "majpower")]
        dyad <- D51[index, "dyad"]
        d <- na.omit(data.frame(Y,X,Z,dyad))
        if (length(unique(X))==1){
            res[b,1] <- NA; res[b,2:3] <- c(NA,NA)
            warning("All obs assigned to one treatment condition at iteration ", b)
        } else {
            fit <- geeglm(Y ~ X + smldep + lcaprat2 + allies + contigkb + logdstab + majpower, 
                family=binomial(link="logit"), data=d, id=dyad, corstr="ar1")
            nd1 <- nd0 <- data.frame(model.matrix(fit)[,-1])
            nd1[,1] <- 1
            nd0[,1] <- 0
            P1 <- predict(fit, nd1, type="response")
            P0 <- predict(fit, nd0, type="response")
            res[b,1] <- mean(pmax(0, P1-P0)/P1, na.rm=T)
            res[b,2:3] <- c(mean(pmax(0, P1-P0)/P1, na.rm=T), 
                mean(pmin(P1, 1-P0)/P1, na.rm=T))
        }
    cat("Threshold ", j, " Resample ", b, " Done. \n")
    }

    res51.post[[j]] <- na.omit(res)

    # Point estimates
    fit.post51 <- geeglm(dispute ~ auth + smldep + lcaprat2 + allies + 
        contigkb + logdstab + majpower, family=binomial(link="logit"), data=D51,
        id=dyad, corstr="ar1")
    nd1.post51 <- nd0.post51 <- data.frame(model.matrix(fit.post51)[,-1])
    nd1.post51[,1] <- 1
    nd0.post51[,1] <- 0
    P1.post51 <- predict(fit.post51, nd1.post51, type="response")
    P0.post51 <- predict(fit.post51, nd0.post51, type="response")

    res51.post.pearl[j] <- mean(pmax(0, P1.post51-P0.post51)/P1.post51, na.rm=T)

    res51.post.exog[j,] <- c(mean(pmax(0, P1.post51-P0.post51)/P1.post51, na.rm=T), 
        mean(pmin(P1.post51, 1-P0.post51)/P1.post51, na.rm=T))

    # Confidence intervals
    CI51.post.pearl[j,] <- c(res51.post.pearl[j] - qnorm(.975)*sd(res51.post[[j]][,1]),
        res51.post.pearl[j] + qnorm(.975)*sd(res51.post[[j]][,1]))

    CI51.post.exog[j,] <- IMCI(l=res51.post.exog[j,1], u=res51.post.exog[j,2],
        var.l=var(res51.post[[j]][,2]), var.u=var(res51.post[[j]][,3]))$ci
}

# Plotting
pdf("fig2.pdf", width=8.5, height=8.5)
par(mfrow=c(2,1), oma=c(1,1,1,1), mar=c(4.1,4.1,2.1,2.1))

point <- res.post.pearl
bounds <- res.post.exog
CI.point <- CI.post.pearl
CI.bounds <- CI.post.exog
plot(-.1,point[1], pch=16, ylim=c(0.3,1), xlim=c(0,10), cex=1.2,
    ylab = "Probability of Causal Attribution", xlab = "",
    main = "Interstate Militarized Disputes, 1886-1992")
lines(c(-.1,-.1), CI.point[1,])
points(c(-.1,-.1), CI.point[1,], pch="-", cex=1.5)
lines(c(.1,.1), bounds[1,], lwd=4)
lines(c(.1,.1), CI.bounds[1,])
points(c(.1,.1), CI.bounds[1,], pch="-", cex=1.5)
for (i in 2:11){
    points(i-1.1, point[i], pch=16, cex=1.2)
    lines(c(i-1.1,i-1.1), CI.point[i,])
    points(c(i-1.1,i-1.1), CI.point[i,], pch="-", cex=1.5)
    lines(c(i-.9,i-.9), bounds[i,], lwd=4)
    lines(c(i-.9,i-.9), CI.bounds[i,])
    points(c(i-.9,i-.9), CI.bounds[i,], pch="-", cex=1.5)
}
abline(h=0.4, lty="dotted")
abline(h=0.6, lty="dotted")
abline(h=0.8, lty="dotted")
abline(h=1, lty="dotted")

point <- res51.post.pearl
bounds <- res51.post.exog
CI.point <- CI51.post.pearl
CI.bounds <- CI51.post.exog
plot(-.1,point[1], pch=16, ylim=c(0.3,1), xlim=c(0,10), cex=1.2,
    ylab = "Probability of Causal Attribution", xlab = "",
    main = "Interstate Militarized Disputes, 1951-1992")
lines(c(-.1,-.1), CI.point[1,])
points(c(-.1,-.1), CI.point[1,], pch="-", cex=1.5)
lines(c(.1,.1), bounds[1,], lwd=4)
lines(c(.1,.1), CI.bounds[1,])
points(c(.1,.1), CI.bounds[1,], pch="-", cex=1.5)
for (i in 2:11){
    points(i-1.1, point[i], pch=16, cex=1.2)
    lines(c(i-1.1,i-1.1), CI.point[i,])
    points(c(i-1.1,i-1.1), CI.point[i,], pch="-", cex=1.5)
    lines(c(i-.9,i-.9), bounds[i,], lwd=4)
    lines(c(i-.9,i-.9), CI.bounds[i,])
    points(c(i-.9,i-.9), CI.bounds[i,], pch="-", cex=1.5)
}
abline(h=0.4, lty="dotted")
abline(h=0.6, lty="dotted")
abline(h=0.8, lty="dotted")
abline(h=1, lty="dotted")
title(sub="Democracy Threshold (Polity Score)", line=3)

dev.off()


